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Abstract 

Two phase Hows can he found in broad situations in nature, bioiogv. and industry devices 
and can involve diverse and complex mechanisms. While the phvsical models may be specific 
for certain dniarions. the mathematical humiliation tnd numerical *:>\itment tor solving die 
governing equations cun be general. Based on the continuum mechanics, we shall treat the fluid 
as a mixrure consisting of two interacting phases -or materials) occupying tin: same region in 
space at any given moment. Hence, we will require information concerning each individual phase 
as needed in a single phase, but also the interactions between them. These interaction terms, 
however, pose additional numerical challenges because they are beyond the basis that we use to 
construct modern numerical schemes, namely the hyperbolicity of equations. Moreover, due to 
disparate differences in time scales, fluid compressibility and nonlinearity become acute, further 
complicating the numerical procedures. In this paper, we will show the ideas and procedure 
how the AUSM-famiiy schemes are extended for solving two phase flows problems. Specifically, 
both phases are assumed in thermodynamic equilibrium, namely, the time scales involved in 
phase interactions are extremely short in comparison with those in fluid speeds and pressure 
fluctuations. Details of the numerical formulation and issues involved are discussed and the 
effectiveness of the method are demonstrated for several industrial examples. 

1 Extension to Real Fluids and Flows with Equilibrium 
Phase Change 

Prior efforts in the construction of AUSM-type algorithms have assumed that the fluid behaves 
as an ideal gas or a mixture thereof. This section details recent extensions of AUSM-type 
schemes that are valid for generalized state equations, which may describe single-phase liquid, 
gas, or supercritical fluid behavior of a given substance. Due to the dramatic differences in 
compressibility among fluids in the different states and possible large differences in the flow 
speed, the “preconditioned” forms of the flux-splitting methods are utilized in the extensions. 
A second thrust of this section is to provide an initial direction toward the development of 
extensions suitable for solving the general multiphase flow problem for arbitrary flow speeds and 
arbitrary levels of compressibility. This initial step starts with the development of an equilibrium 
model for liquid-vapor phase transitions using information extracted from a generalized state 
equation. The resulting equations are similar to the preconditioned, perfect-gas Euler system 
in structure and in mathematical character but may admit such multiphase flow features as 
cavitation zones and vapor-liquid condensation shocks. 

1.1 Real Fluid State Description 

1.1.1 Single Phase Formulation 

The thermodynamic state of a single-phase real fluid is defined by the relations p = p(p, T) and 
h = h(f),T). In the present work, we utilize the Peng-Robinson equation of stare [l], a cubic 
formulation similar to the Van tier Waals equation but generally muen more accurate in the 
liquid phase. The Peng-Robinson equation is given by 

p = Z{p/D,)RT. (IT) 
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Giv**n the fluid pressure and temperature. this equation may he solved for the compressibility 

factor, and the fluid molar density may be then be determined by 


The constant 6 is a function of critical-point fluid properties, while the function a(T) is a function 
of critical point properties and the fluid temperature. The functional form for a(T) can differ 
depending on the polarity of the molecule. The calculations presented herein (for carbon dioxide 
(CO 2 ) and octane) utilize expressions derived in [1] and given below: 


a(T) = a(T c )a(TV,u/), T r = - 

p2'j 1 2 

a(T c ) = 0.45724—^ 

* : 

a(T ri uj) = (1 + *(1 - v^)) 2 

K = 0.37464 + 1.54226a; 

- 0.26992a.' 2 

b = 0.07780^ 


( 1 . 10 ) 

( 1 . 11 ) 


The critical point constants T c and P c and the accentricity factor w are tabulated below for 
CO 2 and octane. 


fluid 

T c (K) 

Pc (Pa) 

CJ 

co 2 

302.2 

73.75e5 

0.225 

octane 

569.4 

24.96e5 

0.4 


The expressions presented above are strictly valid for non-polar molecules. Polar molecules, 
such as water, require different forms for a and 6 to match experimental liquid-state density and 
vapor pressure data properly. 

Typical isotherms for the Peng- Robinson equation are plotted in Fig. 1.1 on a pressure- 
density diagram. Clearly indicated is the vapor regime, where pressure varies nearly linearly 
with density, and the liquid regime, where large pressure changes are required to induce a density 
change. For a given pressure and temperature, the solution of (1.3) returns one or three values of 
the compressibility factor Z, the former corresponding to the single-phase region (either liquid 
or vapor) and the latter corresponding to the two-phase region, where vapor and liquid may 
exist simultaneously. The corresponding densities for a pressure within the two-phase region are 
shown iLs points A-C. A ami C represent saturated vapor and liquid states, while B is physically 
meaningless. For a particular temperature, the '‘allowable” two-phase region is bounded by 
the pressure values at D and E, which are local extrema. The loci of these pressure values 
tor temperatures between the triple and critical points define liquid and vapor spinodal curves. 
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Figure 1.1: Pressure vs. molar density (isotherm below critical temperature). 


dividing the two-phase region into metastable vapor, unstable, and metastable liquid regions. At 
a particular pressure between the liquid and vapor spinodal points, the system is in equilibrium, 
with the vapor and liquid fugacities attaining equal values. This pressure is known as the vapor 
pressure pvap and is calculated as a function of temperature by iterating on the equation 

f(Z v ,T,p) = f(Z t ,T,p), (1.12) 

where the fugacity / is given by 


In — = Z — 1 — ln(Z - B) - 
P 


A Z + (l + y/2)B 
2 V2B V Z + (l->/2)£' 


(1.13) 


The spinodal pressure and density values (points D and E in Fig. LI) can be obtained ana- 
lytically by solving the quartic equation |?|r = 0 , discarding two meaningless roots that occur 
outside the range of validity of the Peng-Robinson equation. The spinodal pressure values bound 
the actual vapor pressure, and an appropriate linear combination can be used as an initial guess 
for the iteration described above. 

The thermodynamic state description for the single-phase fluid is completed by the spec- 
ification of enthalpy departure functions [ 2 ], which introduce a density dependence into the 
enthalpy description and thus account for latent heat effects. For the Peng-Robinson equation, 
the enthalpy per unit mass of the fluid may be expressed as 


h(p, T) = h f (T) + - )-{RT(Z - 1 ) + - ~fjr a(T) ln ( £±il + V ^ )g )|. 

Mw 2\/2b v 7 , /S'n'J 


Z + (l - V2 )B‘ 


(1.14) 


where h f (T) is the enthalpy per unit mass of an ideal gas at the same temperature (determined 
from curve fits presented in McBride, et al [3)). 

The physical sound speed of the fluid can be calculated from thermodynamic considerations; 
a more useful CFD analogue is the acoustic eigenvalue, which may be obtained by determining 
the eigenvalues of the .Jacobian matrix jffc, where F is an Euler flux vector closed according 
to the general expressions p = p(p , T) and ph = ph(p, T). U is the vector of conserved variables, 
and VV is the vector of primitive variables [p, u. u. w. T\ r . The choice of ph (rather than h) and 
the choices of density and temperature as independent variables are dictated by the equilibrium 
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L.L.2 Homogeneous Equilibrium Two-Phase Flow Model 

The Pi*ng- Robinson ’'(juanon (and similar ones! skives no useful information in rfio unsrahle parts 
ol the two-phase region. For densities between the spinodal values, it ran be shown that the 
aroustir eigenvalues are c omplex, meaning that rhe Euler system is not hyperbolic m time and 
that conventional time-marching procedures for integrating the equations are ill posed. It is also 
of note that the liquid spinodal pressure may be negative for high molecular-weight liquids at 
lower temperatures, implying that the simulated expansion of a liquid might produce reasonable 
densities, but unphysical pressures, in the metastable region. 

One means of avoiding cnese difficulties starts with the introduction of a void-fraction formal- 
ism for the two-phase region and the assumption of thermodynamic and kinematic equilibrium 
between the phases. For an equilibrium two-phase flow, the vapor pressure pvap(T) is directly 
related to the temperature through the Clausius-Clapyron equation, and the density and tem- 
perature are independent variables. Given updated values for the density and temperature at a 
grid point as determined from a time-integration method, the following procedure is performed: 

T Determine the vapor pressure at that temperature, either through reference to a curve-fit 
or by the iterative procedure described above, and establish the saturation densities p t (T) 
and p v (T) and the saturation enthalpies h t (T) and k v (T) using (1.3), (1.6) , and (1.14). 

2. If the fluid density is within the saturation limits, the equilibrium equation of state for the 
homogeneous mixture of liquid and vapor is given by 


P — Pvap(T) 

ph(p,T) = p v (T)a v (p 1 T)h v (T) + pt(T)ai(p t T)ht(T) 

Q v(p, T) = f~ p ' (r L 
Pv(T) - pi (T) 

ai(p, T) = 1 - a v 


(1.16) 

(1.17) 

(1.18) 
(1.19) 


3. If the density is not between the saturation values or the temperature is greater than the 
critical value, then the single phase description given by the Peng-Robinson equation will 
be used to determine the pressure and enthalpy. 

In this description, the saturation-state values are strict functions of temperature; density de- 
pendence is introduced through the void fractions a and latent-heat effects arise through the 
change in departure enthalpy between the saturation states. The thermodynamic derivatives 
Pi>’ Pt< (M)p. and (ph) T needed in the time-integration method and in the sound speed defi- 
nition can be computed by straightforward differentiation of the expressions above. These are 
discontinuous at phase transition points, leading to dramatic changes in the effective ‘sound 
speed” in the two-phase region. Figure 1.2 plots a 2 as a function of molar density for both 
the Peng-Robinson equation and the Peng-Robinson equation augmented by the equilibrium 
two-phase flow model (1.16-1.19). The fluid is octane at a temperature of 350 K. As shown, the 
equilibrium two-phase description preserves a real value for the “sound speed”, while the basic 
Peng-Robinson equation results in negative values for a 2 . Also shown is a theoretical result for 
the sound speed in a homogeneous two-phase mixture of liquid and vapor [4]: 

J_ Of 

Pn- ~ p,.n;[T) ~ m j(T ! (I - >,)) 

where nf, t (T ) are obtained from (1.15) evaluated at. the saturation states fj„ t(T). The eigenvalue 
calculation for <r agrees reasonably well with the theoretical estimate except near single-phase 


1 1 




Figure 1.2: a 2 vs. molar density (octane at 350 K). 


/ two-phase junctures, where the latter blends smoothly with the saturation-state values and 
the former exhibits a jump discontinuity. The theoretical expression for the sound speed is 
numerically more robust and is used in all calculations presented herein. Both expressions 
result in very small values (on the order of meters per second) for the “sound speed” near the 
liquid phase / two-phase interface, meaning that a shift to a locally “supersonic” flow condition 
during a phase transition is a distinct possibility. 

The above formulation neglects velocity-slip effects, with the velocity actually solved for 
being a phase-weighted average velocity. This system is hyperbolic in character and is similar 
to the Euler system in structure but admits such multiphase features as cavitation zones and 
condensation shocks. A key element is the use of density and temperature as the “working” 
thermodynamic variables, particularly in contrast with the low-speed formulation described 
earlier, which utilizes pressure and temperature as the “working” thermodynamic variables. 
This choice is driven by the equilibrium closure for the two-phase region, in which pressure and 
temperature are not independent variables. 


1.2 Time-Derivative Preconditioning 

The utility of time-derivative preconditioning in the solution of the real fluid system described 
above lies in its ability to provide a smooth transition between nearly incompressible conditions 
(such as liquid phase or low Mach number vapor or supercritical fluid phase flows) and strongly 
compressible conditions (such as two-phase flows or high-speed vapor phase flows). As discussed 
previously, modifications to AUSM-type discretization are required to extend their range of 
applicability to flows at all speeds. These modifications depend on the choice of preconditioner, 
through the use of the eigenvalues of the preconditioned system. As in the previous work, the 
real fluid extension utilizes the preconditioner of Weiss and Smith [5], which may be expressed 
as a rank-one perturbation of the Jacobian matrix The time-derivative term in the real 
fluid Euler system is replaced by 
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A/; = man L m«ui. A/ ‘ . A/,:„ ; (1.25) 

'.vhnv’ \{ is t\\r lucal Mach number (based on the velocity magnitude i ami A/, : „ is a user-specified 
< urnif Mm h number discussed earlier. Again. the vector IT is ( husen as ; p.n.r.u\ f\ r so that 
die proper closure of the system for an equilibrium two-phase How is achieved. At low speeds, 
this preconditioner essentially replaces the physical thermodynamic derivative p p with 1/U-. 
r°M ■dim' the eigenvalues of the Euler system so that the* eumlitiun number remains bounded. 
This allows uniform convergence in both low-speed and high-speed Hows. The eigenvalues of 
r~‘ A A = OFfdW ) are u .11 ± a . where n is the velocity component hi the x direction and 
u ± a' are as obtained in (??). Note that the real fluid state description does not affect the 
form of the eigenvalues - only the “sound speed” must be redefined. 


1.3 AUSM-type Algorithms for Real Fluids 

Procedures for extending AUSM-type algorithms to operate effectively in conjunction with time- 
derivative preconditioning have been proposed in [6]. These methods reduce to a standard 
upwind formulation at sonic transitions, preserving the discontinuity-capturing traits of the 
methods, but recover viable discretization of the incompressible flow equations as the Mach 
number approaches zero. As the real-fluid state description shares a structured similarity with 
the Euler system with and without preconditioning, it is anticipated that modifications to 
AUSM-type algorithms to allow accurate capturing of real fluid phenomena at all speeds should 
be relatively straightforward. 

A key element in the construction of “all-speed” AUSM-type flux-splitting schemes is the 
need for including a pressure-diffusion term to couple the pressure and velocity fields at low 
Mach number. At higher speeds, the effect of the pressure-diffusion contribution is reduced 
(for AUSMDV and LDFSS) or eliminated (for AUSM-h). For “preconditioned” AUSM-f, the 
pressure-diffusion contribution to the interface flux can be written as 


1/2 


pd, AUSM+ 
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i/2 ausm+ 


Pi ~~ Pi+l 
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P * Pi + 1 
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(1.26) 


where a is the “preconditioned” sound speed defined in (??) and evaluated using averages of left 
(i) and right (1 + 1 state information. The quantity M[ •> ^ ^ interface Mach number 
function, defined as 


^ 1/2 AUSM+ ~ /V, ( + 1)( ; V‘) “ - vt |4.j)(- V ^ J -i) +- v, (i ) (M+l) (1.27) 

The subscript notation i/i + 1 on the vector of ‘advected" variables [1, u, v , w y H\ r indicates its 
evaluation at either the left or the right state, depending on the sign of the complete interface 
mass flux (advective contribution plus pressure diffusion contribution). [6] The^preconditioned” 
version of LDFSS [7] contains a similar term: 
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order unity for gases, but for liquids and supercritical fluids governed bv the Peng-R.obinson 
equation, the coefficient may become much larger than unity. This represents an unphvsical 
source of numerical diffusion for liquid-state calculations, one easily eliminated by redefining 
the pressure-diffusion contributions as 


Fl/2 pd, AUSM+ a ' /2 bv/. 2 1 )' V 'C/2 auSM+ X a? +aj£i 


i/i+1 
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(1.31) 


(1.32) 


To enable an exact reduction to the appropriate incompressible limiting form, the interface 
sound speed can be redefined as 


ai/2 AUSM+ 


-vR 


+ «Li) 


and as 


° 1/2 ldfss - V 


(1.33) 


(1.34) 


All other aspects of the flux-splitting are the same as outlined in Sections ?? and ?? and Refs. 
[6] and [7). It should be noted that the modifications do affect the response of the schemes for 
perfect-gas calculations. For both schemes, the magnitude of the pressure diffusion contribution 
is lowered by a factor of l/y, while for AUSM +, the ability of the scheme to capture a stationary- 
shock wave with no intermediate point is disrupted by the definition in (1.33). Neither of these 
differences affects the performance of the schemes strongly. 

1-4 Higher Order Extension 

To extend the methods outlined above to second-order spatial accuracy, we utilize slope-limited 
Fromm interpolations of the primitive- variable vector [p,u.u,w.T\ r to the t + 1/2 interface. 
As the state description is quite complex and expensive to calculate, some simplifications are 
employed. Fiist, only the pressure and enthalpy are determined from the interpolated density 
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Tn these, the notation L,R represents the use of interpolated values, whereas the notation i, i + 1 
represents the use of nodal values. Results shown later utilize both the minmod and Van Leer 
limiters in defining the interpolated data. 


1.5 Applications 

The techniques outlined in earlier sections have been incorporated into an implicit Navier- 
Stokes solver based on Gauss-Seidel relaxation [8], recently extended to multiblock domains. To 
minimize modifications to the code and to facilitate possible changes in the state description 
the thermodynamic derivatives p„ pr, (ph) p and (ph) T are computed and stored as arrays! 
then used as needed in the construction of the flux Jacobian matrices. Viscosity and thermal 
conductivity data for the single-phase regions are taken from [9] and [10]. For the two-phase 
region, it is assumed that the mixture viscosity and thermal conductivity can be expressed as 
void-fraction weighted averages of the saturation-state values. A quasi- 1-D Euler solver has 
also been written to test basic attributes of the methods. The test cases below illustrate some 
general features of the schemes. 


1.5.1 “Faucet” Problem 

The “faucet” problem [11] is a classic test case for two-fluid codes. In the present context, the 
fluid is taken to be liquid octane in kinematic and thermodynamic equilibrium with its vapor 
at a temperature of 350 K and a vapor pressure of 2061 Pa. The calculation encompasses 
a one-dimensional domain of 12 m, with the inflow conditions specified by the temperature, 
the void fraction of octane vapor (taken as 0.2), and the velocity of the stream (taken as 10 
m/s). The solution is forced by a gravity vector aligned in the direction of the flow, leading to 
acceleration of the fluid and an increase in the vapor-phase void fraction as the density decreases. 
A steady solution is obtained over time, with the transient response being the propagation of 
a discontinuous void wave downstream. With the present closure for the two-phase region, the 
problem i.s hyperbolic in the flow direction, as the effective “sound speed” is much smaller than 
the 10 m/s velocity. As such, all variables are fixed at the inflow and all are extrapolated at 
the outflow. Figure 1.3 presents calculation results for the void wave profile at a particular 
instance m time with ,m analytic solution for the two-phase, separated How problem. A simple 
Euler explicit integration method is used. Only results from AUSM+ are shown, as both schemes 
revert to the same upwind-biased discretization for this problem. The first-order upwind scheme, 
operating at a L'FL of I I), captures the void discontinuity rather sharply but diffuses the peak 



Figure 1.3: Octane vapor void fraction vs. x and time. 


CFT nfrn Second ' order mmm °d and Van Leer - limited extensions require a much smaller 
vf r , ma ! nta,n a reasonat >le de Sree of monotonicity in the void fraction profile. The 
dl “ ' ‘™ ted C3Se P™! ldes a sli « ht improvement in the resolution of the peak value but 
due to the reduction in CFL number, the overall results are only marginally better than the 

a consp der mm “° d - 1,m,ted resu,t is sli g h *iy worse than the first-order result, again 

a consequence of the lower CFL number. It is likely that the second-order results would improve 
witn the use of a more appropriate integration scheme. 

1.5.2 Quasi - 1-D Liquid Expansion 

" a5e L C0nS ‘ derS the flow of initia lly liquid octane through a converging-divemng 
nozzle defined by the area relationship 08 0 g 

A{x) — 1 + 4(i — 1/2)", 0 < x < 1 27) 

The initial conditions are p = 4 x 10 7 Pa T = 340 K and u - in m /e „u . 

n *7 . . , ’ and u - 10 m/s, with the nozzle exit 

thtTZ °' 7 . 4 he imtial P ressure leveI - This problem mimics a cavitating flow in 

enoLh e m P r SUre , drOP eXPenenCed M the fl0W accelerates thr °ugh the nozzle throat is steep 

back in ihl r a ; r r tIOn t0 , the VaP ° r Ph3Se - ThC Rxed Cxit Pressure forces a recompression 
back to the liquid state, sum, fating the collapse of a cavitation region. Figure 1.4 presents 

IZ S Z ? St ?K Utl ° nS f ° r ‘ hree o state equations: the Peng-Robinson (P-R) state description 
pin* r k : qUlhbnum t wo - ph as e flow model, an ideal gas equation of state, and the unmodified 
, f\ ° ins ° n st ate description. As shown, the pressure level in the throat lowers to unphvsical 
° r the ua m° dlfied Peng-Robinson equation, representing a progression into the metastable 
q region- The equation system remains hyperbolic, however, as the density does not drop 
hq “ ld s P lnoda J value. In contrast, the pressure level for the Peng-Robinson equation 
with the equilibrium two-phase flow model lowers to the vapor pressure of octane. This results 
in the generation of a vapor phase and a decrease in the fluid density (Fig. 1.5). As expected 
e 'quid octane density varies little in the convergent section of the nozzle. The shock-like 
ASr+TnTr rwo-phase fluid back to the liquid state is captured well by both the 

. .' . SS discretizations. Some effects of the higher-order extensions can be seen in 

° P os ' t,0, “ug »f the condensation shock, and little difference between the AUSM+ and LDFSS 
pm ictions is observed. Pressure distributions from a perfect-gas closure are also shown in Fig. 
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Figure 1.4: Pressure vs. x: expanding liquid octane. 


th * COmp Z' S T ^ ighIi « htin « the «Pected differences in the flow response due to the 
state description. The shock wave is captured monotonically by all methods. 



Figure 1.5: Density and vapor void fraction vs. x: expanding liquid octane. 


1.5.3 


Liquid CO, Expansion Through a Sharp Orifice Nozzle 


gure 1.6 illustrates axial velocity and density contours in the interior of a reservoir / capillary 
nozzle system for spraying liquid CO , A two- block grid is utilized, with the reservoir block 
containing 6oxto3 points and the capillary nozzle block containing 97x97 points The flow is 
axisymmetnc, and the reservoir total conditions ,ire p„ = 10 x I0 a Pa and T 0 = 290A' These 
conditions place the incoming fluid in the liquid state. The inflow boundary conditions fix the 
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Figure 1.6: Velocity and density contours: liquid CO 2 expansion (AUSM+ upwinding). 

No discontinuity in the bulk fluid properties is evidenced, however, as the transition into the 
two-phase region takes place very near the critical point of the fluid. Calculations without the 
equilibrium two-phase flow model were found to be unstable, as the rapid expansion drives the 
fluid density in the corner region below the liquid spinodal value. The pressure values remain 
reasonable, in contrast to the octane expansion described above, but the acoustic speeds become 
complex. Only AUSM+ solutions are shown; LDFSS solutions are very similar. The effects of 
the mmmod-limited second-order extension are confined to the orifice region, where the second- 
order calculation results in more crisp predictions of the supersonic flow response. Features 
of note include regular oblique-shock and Mach wave reflections as well as a small pocket of 
reversed flow downstream of the corner. The displacement effect of this structure forms an area 
minimum, allowing the transition to supersonic flow. 


1.5.4 Liquid Octane Expansion Through a Sharp Orifice Nozzle 

The fourth test case involves the acceleration of liquid octane through a capillarv tube. Devices 
similar to this are used in fuel injection systems. Again, a two-block grid is considered, with 
reservoir containing 65x153 nodes and the capillary tube containing 97x97 nodes. This problem 
is also axisymmetric, and the reservoir conditions are p 0 = 10 x 10 6 Pa and T 0 = -400 AT. Octane is 
liquid under these conditions. The closeup in Fig. 1.7 plots density contours in the vicinity of the 
reservoir / tube juncture. The rapid pressure drop experienced as the fluid accelerates around 
the sharp corner cavitates the fluid, initially producing a bubble of nearly pure vapor. The wake 
if the bubble is a t vv o- p h.iso mixtiiie of fluid, characterized by an increasing liquid content as 
the flow proceeds downstream. The two-phase / single-phase interface is sharplv captured near 
the corner, with the shape of the bubble determined by a balance between the pressure jump 
and numerical surface tension resulting from the upwinding. Differences between first-order 
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Figure 1.7: Density contours: octane cavitation. 


AUSM+ and LDFSS solutions are minimal. These solutions are extremely difficult to obtain. 

ith the current methods, the iteration must be nearly time-accurate with global time-stepping 
in order to capture stable bubble growth. Local time steps result in the intermittent appearance 
of octane vapor bubbles, which grow, propagate, and collapse in a highly unsteady manner With 
the present thermodynamic model, the collapse of a vapor bubble results in higher temperatures 
raising the vapor pressure and promoting more bubble growth. While many of these trends are 
consistent with the physics of cavitation bubble formation [4], the robustness and efficiency of 
the current procedures in capturing steady bubble behavior is a concern. The same is true for 
e higher-order extensions, which are even more susceptible to transient bubble growth As a 
result, no higher-order solutions for this problem are yet available. Cavitation calculations using 
the Sanchez-Lacombe [12J equation of state, a lattice-fluid formulation valid for high molecular 
weight liquids, are underway for octane and water; these indicate somewhat better numerical 


1.6 Concluding Remarks 

fiZd T difiC , at r S f0r e , Xtending AUSM+ and LDFSS low-diffusion upwind schemes schemes 
the calculation of real fluids at all speeds and at all states of compressibility have 
been outlined in this section. The real fluid state description is based on the Peng-Robinson 
equation, enhanced by an equilibrium model for liquid-vapor phase transitions. Results indicate 
that the modifications proposed herein are effective in simulating incompressible liquid and 
mpressible vapor responses as well as multiphase flow phenomena, such the appearance of 
cavitation bubbles and vapor-liquid condensation shocks. A point of concern is the robustness 
ot the current procedures capturing stationary cavitation bubbles - mollifications to improve 
this behavior are underway. This work provides a starting point for a more comprehensive 
investigation of upwind discretization techniques for general nonequilibrium multiphase flows - 
• ‘Itotf.s in this dirrctiuii .in* ,ilso un<ItTw;iv 
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